Considering the prevalence of obesity (1), its associated fracture risk [(2)] and that these patients are especially prone to develop loading associated musculoskeletal injuries [(5)], there is a need to develop mechanical loading prediction equations that are accurate for overweight and obese subjects in order to precisely determine and monitor exercise-associated mechanical loading in these patients. Therefore, the purpose of this study was to develop prediction equations based on accelerometer data able to accurately predict ground reaction forces (GRF) and loading rate (LR) on a broad range of body masses, from normal weight to severely obese subjects, setting thereby, a base for objective prescription and monitoring of exercise mechanical loads.
library(knitr)
library(broman)
library(tidyverse)
library(here)
library(nlme)
library(piecewiseSEM)
library(ez)
source(here("R", "BMI_categories.R"))
source(here("R", "cross_validate_mixed_model.R"))
source(here("R", "accuracy_indices.R"))
source(here("R", "get_BA_plot.R"))Initially, 4 files needed to be loaded, containing data from each of the accelerometer placements utilized in the study (lower back and hip), and for each of the mechanical loading variables (GRF and LR).
back <- read_csv(here("data", "back_sec.csv")) %>%
select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>%
filter(speed <= 6)
hip <- read_csv(here("data", "hip_sec.csv")) %>%
select(ID, speed, body_mass, height, BMI, sex, age, pVGRF_N, pVACC_g, pRGRF_N, pRACC_g) %>%
filter(speed <= 6)
anthopometric <- read_csv(here("data", "back_sec.csv")) %>%
select(ID, speed, height, BMI, sex, age)
back_LR <- read_csv(here("data", "max_rates_back_IMU.csv")) %>%
left_join(anthopometric, by = c("ID", "speed")) %>%
BMI_categories() %>%
select(
ID, speed, body_mass, height, BMI, BMI_cat, sex, age, pVLR_Ns, pVATR_gs, pRLR_Ns, pRATR_gs
)
hip_LR <- read_csv(here("data", "max_rates_hip_IMU.csv")) %>%
left_join(anthopometric, by = c("ID", "speed")) %>%
BMI_categories() %>%
select(
ID, speed, body_mass, height, BMI, BMI_cat, sex, age, pVLR_Ns, pVATR_gs, pRLR_Ns, pRATR_gs
)A piece of a data frame used in all subsequent analysis can be seen below.
## # A tibble: 300 x 11
## ID speed body_mass height BMI sex age pVGRF_N pVACC_g pRGRF_N
## <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 3 2 85.6 171 29.3 M 50 1023. 1.26 1030.
## 2 3 3 85.6 171 29.3 M 50 1118. 1.51 1133.
## 3 3 4 85.6 171 29.3 M 50 1202. 1.61 1230.
## 4 3 5 85.6 171 29.3 M 50 1200. 1.63 1246.
## 5 3 6 85.6 171 29.3 M 50 1232. 1.73 1294.
## 6 4 2 61.1 162 23.3 F 45 671. 1.16 678.
## 7 4 3 61.1 162 23.3 F 45 702. 1.23 716.
## 8 4 4 61.1 162 23.3 F 45 756. 1.35 770.
## 9 4 5 61.1 162 23.3 F 45 809. 1.50 828.
## 10 4 6 61.1 162 23.3 F 45 912. 1.67 926.
## # … with 290 more rows, and 1 more variable: pRACC_g <dbl>
## # A tibble: 271 x 12
## ID speed body_mass height BMI BMI_cat sex age pVLR_Ns pVATR_gs
## <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> <dbl> <dbl> <dbl>
## 1 3 2 86.5 171 29.3 overwe… M 50 4367. 7.95
## 2 3 3 86.5 171 29.3 overwe… M 50 7895. 13.3
## 3 3 4 86.5 171 29.3 overwe… M 50 10585. 14.1
## 4 3 5 86.5 171 29.3 overwe… M 50 12798. 15.1
## 5 3 6 86.5 171 29.3 overwe… M 50 13556. 20.2
## 6 4 2 61.1 162 23.3 normal… F 45 3650. 4.12
## 7 4 3 61.1 162 23.3 normal… F 45 5049. 5.42
## 8 4 4 61.1 162 23.3 normal… F 45 6872. 11.3
## 9 4 5 61.1 162 23.3 normal… F 45 7992. 17.4
## 10 4 6 61.1 162 23.3 normal… F 45 10800. 15.3
## # … with 261 more rows, and 2 more variables: pRLR_Ns <dbl>,
## # pRATR_gs <dbl>
We first determined the relationship between peak acceleration (pACC) and peak GRF (pGRF), and between peak ACC acceleration transient rate (pATR) and peak LR (pLR) obtained during the incremental walking speeds used in our experimental protocol. This was performed in order to verify the consistency of these relationships for increasing pACC and pATR values and the ability of the accelerometer to discriminate differences in pGRF and pLR between subjects in different body mass indices (BMI) classes. A scatterplot with these relationships for all two accelerometer placements tested and for both the resultant and its vertical component is depicted in Figure 1. Generally, as expected, pACC recorded by accelerometers in all placements, showed a linear increase as a function of pGRF increases. Also, the recorded ACC were shown to be able to provide a good discrimination between different BMI classes, as for the same registered pACC the pGRF tended to be consistently higher for subjects in higher BMI classes. As was also expected, the pATR recorded by the accelerometers demonstrated a linear increase as the pLR increases. Contrarily to what was observed on the relationship between pACC and pGRF, the pATR x pLR plot does not provide a good discrimination between different BMI classes. Moreover, the values at low magnitudes showed a high degree of clustering, whereas the values at higher magnitudes tended to a greater dispersion.
Code to generate these plots can be seen here.
In a first attempt to predict the pGRF using accelerometry data, Newton’s second law of motion (force = mass · acceleration) was applied to predict peak resultant GRF (pRGRF) from peak resultant ACC (pRACC) derived from each accelerometer placement. As both accelerometers were located at waistline, both were close to the body center of mass, and therefore, their recorded ACC values could represent the ACC of the body center of mass.
The code usede to make these predictions can be seen below.
# Back
back_2nd_law <- back %>%
mutate(
pRACC_ms2 = pRACC_g * 9.81,
pRGRF_N_predicted = body_mass * pRACC_ms2
)
# Hip
hip_2nd_law <- hip %>%
mutate(
pRACC_ms2 = pRACC_g * 9.81,
pRGRF_N_predicted = body_mass * pRACC_ms2
)To evaluate these predictions, some accuracy indices were computed. These indices were the mean absolute error (MAE), mean absolute percent error (MAPE) and root mean square error (RMSE), and they were computed with the accuracy_indices() function.
# Back
back_2nd_law_accuracy <- accuracy_indices(back_2nd_law, "pRGRF_N", "pRGRF_N_predicted")
back_2nd_law_accuracy## bias_mean bias_sd LoA_inf LoA_sup MAE MAE_sd RMSE bias_mean_p
## 1 -90.37 137.71 -360.28 179.54 109.86 122.68 164.52 -7.21
## bias_sd_p LoA_inf_p LoA_sup_p MAPE MAPE_sd CV_RMSE
## 1 10.7 -28.18 13.76 9.21 9.03 14.06
# Hip
hip_2nd_law_accuracy <- accuracy_indices(hip_2nd_law, "pRGRF_N", "pRGRF_N_predicted")
hip_2nd_law_accuracy## bias_mean bias_sd LoA_inf LoA_sup MAE MAE_sd RMSE bias_mean_p
## 1 -224.65 191.6 -600.19 150.89 229.4 185.86 295.05 -18.31
## bias_sd_p LoA_inf_p LoA_sup_p MAPE MAPE_sd CV_RMSE
## 1 13.86 -45.48 8.86 18.85 13.11 25.25
Linear mixed models (LMM) were developed to predict peak resultant GRF (pRGRF) and peak vertical GRF (pVGRF). Distinct LMMs were developed with data from lowee back and hip accelerometer placement (back and hip data frames) using the lme() function of the nlme package. Covariance structure used was an autoregressive process of order 1 (correlation = corAR1) and maximum likelihood method was used for estimating parameters (method = "ML"). Predictors tested on pRGRF models were body mass and peak resultant acceleration (pRACC), while on pVGRF were body mass and vertical acceleration (pVACC). All of them were tested as fixed effects and have shown to be significant predictors (e.g.: fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass). Both random intercept and slopes were tested, but only the random intercept inclusion has showed models improvement (random = ~ 1 | ID). Linear, quadratic and cubic polynomial simulations were also tested, whereas the last one did not contribute significantly to the models. Final models were chosen according to -2 log-likelihood statistics. Traditional coefficient of determination (R2) was represented by conditional R2 (6) computed with rsquared() function of the piecewiseSEM package, that estimates the variance explained by the whole model (7).
Code to build the LMMs, as well as their summary and R2, can be found below.
# For resultant peak ground reaction force
back_res_LMM <- lme(
fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back
)
r2_back_res_LMM <- rsquared(back_res_LMM)
summary(back_res_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back
## AIC BIC logLik
## 3185.997 3215.628 -1584.999
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.0184402 79.53301
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8530994
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -698.7031 123.44832 233 -5.659884 0e+00
## pRACC_g 1047.5129 122.40109 233 8.558036 0e+00
## I(pRACC_g^2) -345.2605 36.10073 233 -9.563809 0e+00
## body_mass 3.8294 1.06086 62 3.609729 6e-04
## pRACC_g:body_mass 6.0219 0.61893 233 9.729449 0e+00
## Correlation:
## (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g -0.896
## I(pRACC_g^2) 0.648 -0.892
## body_mass -0.636 0.268 0.143
## pRACC_g:body_mass 0.565 -0.281 -0.167 -0.922
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.0955967 -0.6001009 -0.0930397 0.5792458 3.1083008
##
## Number of Observations: 300
## Number of Groups: 64
## Response family link method Marginal Conditional
## 1 pRGRF_N gaussian identity none 0.9361655 0.9361655
hip_res_LMM <- lme(
fixed = pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip
)
r2_hip_res_LMM <- rsquared(hip_res_LMM)
summary(hip_res_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip
## AIC BIC logLik
## 3091.034 3120.637 -1537.517
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02080893 71.4172
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8662087
## Fixed effects: pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass + pRACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -300.9909 89.71869 232 -3.354829 9e-04
## pRACC_g 522.6850 71.61994 232 7.298037 0e+00
## I(pRACC_g^2) -171.5606 16.91100 232 -10.144910 0e+00
## body_mass 3.9596 0.82234 62 4.814986 0e+00
## pRACC_g:body_mass 5.3671 0.41930 232 12.799992 0e+00
## Correlation:
## (Intr) pRACC_g I(RACC bdy_ms
## pRACC_g -0.878
## I(pRACC_g^2) 0.588 -0.854
## body_mass -0.766 0.422 0.030
## pRACC_g:body_mass 0.678 -0.471 -0.038 -0.891
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.54273847 -0.59285086 0.05857753 0.70507501 2.51256437
##
## Number of Observations: 299
## Number of Groups: 64
## Response family link method Marginal Conditional
## 1 pRGRF_N gaussian identity none 0.9488446 0.9488446
# For vertical peak ground reaction force
back_vert_LMM <- lme(
fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back
)
r2_back_vert_LMM <- rsquared(back_vert_LMM)
summary(back_vert_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back
## AIC BIC logLik
## 3234.998 3264.628 -1609.499
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02508531 102.5201
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.9079104
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -776.8934 145.57197 233 -5.336834 0
## pVACC_g 1042.9052 146.96000 233 7.096524 0
## I(pVACC_g^2) -336.2115 44.82805 233 -7.500024 0
## body_mass 6.2132 1.22304 62 5.080120 0
## pVACC_g:body_mass 5.0805 0.72861 233 6.972936 0
## Correlation:
## (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g -0.889
## I(pVACC_g^2) 0.652 -0.894
## body_mass -0.688 0.327 0.066
## pVACC_g:body_mass 0.587 -0.341 -0.101 -0.891
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.2287736 -0.6177471 -0.1216603 0.4585545 3.7585643
##
## Number of Observations: 300
## Number of Groups: 64
## Response family link method Marginal Conditional
## 1 pVGRF_N gaussian identity none 0.8935621 0.8935621
hip_vert_LMM <- lme(
fixed = pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g : body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip
)
r2_hip_vert_LMM <- rsquared(hip_vert_LMM)
summary(hip_vert_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip
## AIC BIC logLik
## 3153.414 3183.017 -1568.707
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.02202366 84.32584
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.8869442
## Fixed effects: pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass + pVACC_g:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -435.7365 122.69794 232 -3.551294 5e-04
## pVACC_g 586.6627 112.66171 232 5.207295 0e+00
## I(pVACC_g^2) -188.9689 28.55314 232 -6.618146 0e+00
## body_mass 5.8047 0.97485 62 5.954511 0e+00
## pVACC_g:body_mass 4.9544 0.54067 232 9.163472 0e+00
## Correlation:
## (Intr) pVACC_g I(VACC bdy_ms
## pVACC_g -0.910
## I(pVACC_g^2) 0.715 -0.907
## body_mass -0.759 0.471 -0.125
## pVACC_g:body_mass 0.688 -0.524 0.135 -0.888
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.95040403 -0.53289689 -0.04360304 0.45305786 4.09646433
##
## Number of Observations: 299
## Number of Groups: 64
## Response family link method Marginal Conditional
## 1 pVGRF_N gaussian identity none 0.9263472 0.9263473
In all models, conditional R2 values ranged from 0.89 to 0.95.
Model validation was assessed by the leave-one-out cross-validation (LOOCV) method. This approach is recommended when a different sample is not available for cross-validation and it provides an insight on the model potential to predict outcomes in a new independent sample (8). For LOOCV each participant’s data was separated into a testing dataset (one at a time) with the remaining data being in the training dataset. New LMM, with the same outcomes and predictors as determined for the entire sample, were developed using the training dataset and then used to predict pGRF for the only participant in the testing dataset. This process was repeated for all participants (64 times). Data from testing dataset were used in the remaining statistical analysis.
LOOCV of the LMM was done with the cross_validate_mixed_model() function and is shown below.
# For resultant peak ground reaction force
fix_eff <- pRGRF_N ~ pRACC_g + I(pRACC_g^2) + body_mass
rand_eff <- ~ 1 | ID
# Back
LOOCV_back_res_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_res_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))
# For vertical peak ground reaction force
fix_eff <- pVGRF_N ~ pVACC_g + I(pVACC_g^2) + body_mass
rand_eff <- ~ 1 | ID
# Back
LOOCV_back_vert_LMM <- do.call(rbind, (lapply(unique(back$ID), cross_validate_mixed_model, df = back)))
# Hip
LOOCV_hip_vert_LMM <- do.call(rbind, (lapply(unique(hip$ID), cross_validate_mixed_model, df = hip)))Bland-Altman plots were used to examine the agreement between pGRF measured with force plates and those predicted through the regression equations. The difference of the actual and predicted pGRF was plotted against their mean. Bias was expressed as the mean of these differences and the limits of agreement were obtained using ±1.96 standard deviation of the mean between actual and predicted pGRF (9).
Bland-Altman plots can be seen below for all accelerometer placements. Panel A shows plots of pRGRF and panel B shows plots of pVGRF.
source(here("figs", "fig2.R"))
BA_plot_grid_1 <- plot_grid(
back_res_BA_plot + theme(legend.position = "none"),
back_vert_BA_plot + theme(legend.position = "none"),
hip_res_BA_plot + theme(legend.position = "none"),
hip_vert_BA_plot + theme(legend.position = "none"),
labels = c("A", "B", "", ""),
align = "h", vjust = 1, label_size = 16,
ncol = 2, nrow = 2
)
legend <- get_legend(back_res_BA_plot)
BA_plot_grid <- plot_grid(BA_plot_grid_1, legend, ncol = 1, rel_heights = c(1, 0.1))
BA_plot_gridCode to generate these plots can be seen here.
One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.
# For resultant peak ground reaction force
# Back
LOOCV_back_res_LMM$diff <- LOOCV_back_res_LMM$pRGRF_N - LOOCV_back_res_LMM$pRGRF_N_predicted
LOOCV_back_res_LMM$mean <- (LOOCV_back_res_LMM$pRGRF_N + LOOCV_back_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_back_res_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_back_res_LMM$diff
## t = 0.44993, df = 299, p-value = 0.6531
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.724802 12.303988
## sample estimates:
## mean of x
## 2.289593
# Hip
LOOCV_hip_res_LMM$diff <- LOOCV_hip_res_LMM$pRGRF_N - LOOCV_hip_res_LMM$pRGRF_N_predicted
LOOCV_hip_res_LMM$mean <- (LOOCV_hip_res_LMM$pRGRF_N + LOOCV_hip_res_LMM$pRGRF_N_predicted) / 2
t.test(LOOCV_hip_res_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_hip_res_LMM$diff
## t = 0.36236, df = 298, p-value = 0.7173
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7.774785 11.284133
## sample estimates:
## mean of x
## 1.754674
# For vertical peak ground reaction force
# Back
LOOCV_back_vert_LMM$diff <- LOOCV_back_vert_LMM$pVGRF_N - LOOCV_back_vert_LMM$pVGRF_N_predicted
LOOCV_back_vert_LMM$mean <- (LOOCV_back_vert_LMM$pVGRF_N + LOOCV_back_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_back_vert_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_back_vert_LMM$diff
## t = -0.10563, df = 299, p-value = 0.9159
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -12.19819 10.95537
## sample estimates:
## mean of x
## -0.6214066
# Hip
LOOCV_hip_vert_LMM$diff <- LOOCV_hip_vert_LMM$pVGRF_N - LOOCV_hip_vert_LMM$pVGRF_N_predicted
LOOCV_hip_vert_LMM$mean <- (LOOCV_hip_vert_LMM$pVGRF_N + LOOCV_hip_vert_LMM$pVGRF_N_predicted) / 2
t.test(LOOCV_hip_vert_LMM$diff, mu = 0)##
## One Sample t-test
##
## data: LOOCV_hip_vert_LMM$diff
## t = 0.0087888, df = 298, p-value = 0.993
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -9.956701 10.046032
## sample estimates:
## mean of x
## 0.04466578
Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.
Results showed a significant proportional bias (p < 0.05), however, with a low magnitude (highest R2 = 0.032). These results suggest that despite there is a trend for underestimation at increasingly higher pGRF values, the magnitude of this effect is neglectable.
# For resultant peak ground reaction force
# Back
back_res_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_res_LMM)
summary(back_res_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_res_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -267.995 -55.071 -2.505 52.450 294.243
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -47.77144 19.69289 -2.426 0.01587 *
## mean 0.04282 0.01628 2.630 0.00899 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.28 on 298 degrees of freedom
## Multiple R-squared: 0.02268, Adjusted R-squared: 0.0194
## F-statistic: 6.915 on 1 and 298 DF, p-value: 0.008991
##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_res_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -268.721 -53.909 5.937 51.546 235.762
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -40.38734 18.75117 -2.154 0.0321 *
## mean 0.03610 0.01552 2.325 0.0207 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 83.12 on 297 degrees of freedom
## Multiple R-squared: 0.01788, Adjusted R-squared: 0.01457
## F-statistic: 5.406 on 1 and 297 DF, p-value: 0.02074
# For vertical peak ground reaction force
# Back
back_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_vert_LMM)
summary(back_vert_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_back_vert_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -342.85 -64.14 -4.92 51.55 354.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -46.47386 23.16029 -2.007 0.0457 *
## mean 0.03998 0.01954 2.046 0.0416 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 101.4 on 298 degrees of freedom
## Multiple R-squared: 0.01386, Adjusted R-squared: 0.01055
## F-statistic: 4.187 on 1 and 298 DF, p-value: 0.04162
# Hip
hip_vert_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_vert_LMM)
summary(hip_vert_BA_plot_LR)##
## Call:
## lm(formula = diff ~ mean, data = LOOCV_hip_vert_LMM)
##
## Residuals:
## Min 1Q Median 3Q Max
## -255.386 -53.520 -0.846 49.427 295.833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -42.10744 19.93510 -2.112 0.0355 *
## mean 0.03683 0.01685 2.186 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.33 on 297 degrees of freedom
## Multiple R-squared: 0.01583, Adjusted R-squared: 0.01252
## F-statistic: 4.778 on 1 and 297 DF, p-value: 0.02961
To evaluate models prediction accuracy, MAE, MAPE and RMSE were computed. These indices were computed with the accuracy_indices() function.
# For resultant peak ground reaction force
# Back
back_res_accuracy <- accuracy_indices(LOOCV_back_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# Hip
hip_res_accuracy <- accuracy_indices(LOOCV_hip_res_LMM, "pRGRF_N", "pRGRF_N_predicted")
# For vertical peak ground reaction force
# Back
back_vert_accuracy <- accuracy_indices(LOOCV_back_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")
# Hip
hip_vert_accuracy <- accuracy_indices(LOOCV_hip_vert_LMM, "pVGRF_N", "pVGRF_N_predicted")Table below shows these indices for each accelerometer placement in both pRGRF and pVGRF.
accuracy <- data.frame(
acc_placement = c(rep("Back", 2), rep("Hip", 2)),
condition = rep(c("pRGRF", "pVGRF"), 2),
MAE = c(
round(back_res_accuracy[[5]], 1), round(back_vert_accuracy[[5]], 1),
round(hip_res_accuracy[[5]], 1), round(hip_vert_accuracy[[5]], 1)
),
MAPE = c(
paste(round(back_res_accuracy[[12]], 1), "%", sep = ""), paste(round(back_vert_accuracy[[12]], 1), "%", sep = ""),
paste(round(hip_res_accuracy[[12]], 1), "%", sep = ""), paste(round(hip_vert_accuracy[[12]], 1), "%", sep = "")
),
RMSE = c(
round(back_res_accuracy[[7]], 1), round(back_vert_accuracy[[7]], 1),
round(hip_res_accuracy[[7]], 1), round(hip_vert_accuracy[[7]], 1)
)
)
kable(accuracy, col.names = c("Accelerometer Placement", "GRF", "MAE", "MAPE", "RMSE"), align = c("l", "l", "r", "r", "r"))| Accelerometer Placement | GRF | MAE | MAPE | RMSE |
|---|---|---|---|---|
| Back | pRGRF | 66.7 | 5.9% | 88.0 |
| Back | pVGRF | 74.9 | 6.5% | 101.7 |
| Hip | pRGRF | 65.5 | 5.9% | 83.6 |
| Hip | pVGRF | 65.6 | 5.9% | 87.7 |
All conditions showed a MAE from near 65 N to near 81 N, MAPE from near 6% to near 7% and RMSE from near 84 N to 84 N.
A series of repeated measures analysis of variance (ANOVA) were run to assess whether pGRF predicted from the regression equations were significantly different from those measured with FP. Walking speeds, accelerometer placements (lower back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis. These procedures were taken separately for resultant and its vertical component. If assumptions of sphericity were violated (p < 0.05), the conservative Greenhouse–Geisser correction factor would be applied to adjust the degrees of freedom. Post-hoc analyses would be conducted using pairwise comparisons with Holm’s correction if a significant difference were observed among actual and predicted pGRF.
In order to run the repeated measures ANOVA, new data frames needed to be built, for both pRGRF and pVGRF, putting all the pGRF values in a single variable and grouping then by accelerometer placement or actual value in another variable. An example of code to tidy the data is shown below.
## Predicted pRGRF
back_res_pred <- LOOCV_back_res_LMM %>%
select(ID, speed, pRGRF_N_predicted) %>%
spread(key = speed, value = pRGRF_N_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = back
)
hip_res_pred <- LOOCV_hip_res_LMM %>%
select(ID, speed, pRGRF_N_predicted) %>%
spread(key = speed, value = pRGRF_N_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = hip
)
res_pred_df <- back_res_pred %>%
full_join(hip_res_pred, by = c("ID", "speed")) %>%
na.omit()
## Actual pRGRF
back_res_actual <- LOOCV_back_res_LMM %>%
select(ID, speed, pRGRF_N) %>%
spread(key = speed, value = pRGRF_N) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_back
)
hip_res_actual <- LOOCV_hip_res_LMM %>%
select(ID, speed, pRGRF_N) %>%
spread(key = speed, value = pRGRF_N) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_hip
)
res_actual_df <- back_res_actual %>%
full_join(hip_res_actual, by = c("ID", "speed")) %>%
na.omit()
res_actual_df <- res_actual_df %>%
mutate(actual = (actual_back + actual_hip) / 2) %>%
select(ID, speed, actual)
## Merge predicted and actual data frames
res_ANOVA_df <- res_pred_df %>%
full_join(res_actual_df, by = c("ID", "speed")) %>%
gather(
back, hip, actual,
key = "group",
value = "pRGRF"
)
res_ANOVA_df$ID <- as.factor(res_ANOVA_df$ID)
res_ANOVA_df$speed <- as.factor(res_ANOVA_df$speed)
res_ANOVA_df$group <- as.factor(res_ANOVA_df$group)Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Summary statistics for pRGRF and pVGRF ANOVA can be seen below, respectively.
# For resultant peak ground reaction force
res_ANOVA <- ezANOVA(
data = res_ANOVA_df,
dv = pRGRF,
wid = ID,
within = .(speed, group),
detailed = TRUE,
type = 3
)
res_ANOVA## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 45 8.401226e+08 41784505.3 9.047736e+02 1.922348e-31
## 2 speed 4 180 7.469909e+06 465997.2 7.213474e+02 1.330059e-109
## 3 group 2 90 8.151181e+01 876689.8 4.183956e-03 9.958250e-01
## 4 speed:group 8 360 1.474347e+04 331353.9 2.002259e+00 4.530694e-02
## p<.05 ges
## 1 * 9.508154e-01
## 2 * 1.466746e-01
## 3 1.875619e-06
## 4 * 3.391386e-04
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 speed 0.0353185983 1.104234e-26 *
## 3 group 0.4299955589 8.633907e-09 *
## 4 speed:group 0.0005121196 3.532794e-47 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 speed 0.4083488 2.563967e-46 * 0.4216157 9.726767e-48
## 3 group 0.6369409 9.746061e-01 0.6471940 9.758931e-01
## 4 speed:group 0.2655135 1.378260e-01 0.2790290 1.350822e-01
## p[HF]<.05
## 2 *
## 3
## 4
# For vertical peak ground reaction force
vert_ANOVA <- ezANOVA(
data = vert_ANOVA_df,
dv = pVGRF,
wid = ID,
within = .(speed, group),
detailed = TRUE,
type = 3
)
vert_ANOVA## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 (Intercept) 1 45 8.079216e+08 38787197.4 937.33175604 8.997743e-32
## 2 speed 4 180 5.570015e+06 377459.7 664.04624137 1.443126e-106
## 3 group 2 90 1.426599e+03 1353453.5 0.04743195 9.536992e-01
## 4 speed:group 8 360 7.412484e+03 339268.4 0.98317970 4.487240e-01
## p<.05 ges
## 1 * 9.518633e-01
## 2 * 1.199726e-01
## 3 3.491533e-05
## 4 1.813905e-04
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 speed 0.0449465515 1.555112e-24 *
## 3 group 0.8418464231 2.265287e-02 *
## 4 speed:group 0.0005848457 4.244722e-46 *
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF]
## 2 speed 0.4411131 1.738440e-48 * 0.4577530 3.227609e-50
## 3 group 0.8634433 9.346655e-01 0.8947448 9.396637e-01
## 4 speed:group 0.3093066 3.915361e-01 0.3286990 3.951765e-01
## p[HF]<.05
## 2 *
## 3
## 4
For either resultant and its vertical component, repeated measures ANOVA demonstrated that actual and predicted pGRF increased significantly (p < 0.001) along with increments in walking speed.
To assess for any significant differences between actual and predicted pGRF in each speed, five separate repeated measures ANOVA were done, one for each walking speed.
# For resultant peak ground reaction force
# 2 km/h
res_ANOVA_s2 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 2),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 3 km/h
res_ANOVA_s3 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 3),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 4 km/h
res_ANOVA_s4 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 4),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 5 km/h
res_ANOVA_s5 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 5),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 6 km/h
res_ANOVA_s6 <- ezANOVA(
data = res_ANOVA_df %>% filter(speed == 6),
dv = pRGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# For vertical peak ground reaction force
# 2 km/h
vert_ANOVA_s2 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 2),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 3 km/h
vert_ANOVA_s3 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 3),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 4 km/h
vert_ANOVA_s4 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 4),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 5 km/h
vert_ANOVA_s5 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 5),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)
# 6 km/h
vert_ANOVA_s6 <- ezANOVA(
data = vert_ANOVA_df %>% filter(speed == 6),
dv = pVGRF,
wid = ID,
within = group,
detailed = TRUE,
type = 3
)No significant differences (p > 0.05) between actual and predicted pGRF in each speed were found.
These results can be seen in figure below.
Our equation was compared with a previously published reference equation (11), in which Neugebauer et al. used a similar approach for the prediction of pGRF. Comparison was performed using the equation to predict pVGRF from hip-worn accelerometers, as it was the only suitable with our data. The analysis was performed in three ways using: i) a subsample of participants with normal weight or overweight (BMI ≥ 18.5 and < 30 kg.m-2); ii) a subsample of obese participants (BMI ≥ 30 kg.m-2); iii) the whole sample.
New data frames needed to be built for each of the subsamples and the whole sample. Then, pVGRF was predicted using the reference equation. Code to build the data frames is shown below.
# Prepare data frames
non_obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>%
filter(BMI < 30)
obese_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted) %>%
filter(BMI >= 30)
whole_sample_df <- read_csv("~/Dropbox/Projects/walking_GRF_ACC/LOOCV_hip_vert.csv") %>%
select(ID, speed, body_mass, BMI, BMI_cat, pVACC_g, pVGRF_N, pVGRF_N_predicted)
# Apply Neugebauer 2014 equation
non_obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(non_obese_df)) {
non_obese_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * non_obese_df$pVACC_g[i]) + (0.014 * non_obese_df$body_mass[i]))
}
obese_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(obese_df)) {
obese_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * obese_df$pVACC_g[i]) + (0.014 * obese_df$body_mass[i]))
}
whole_sample_df$pVGRF_N_Neugebauer <- NA
for (i in 1:nrow(whole_sample_df)) {
whole_sample_df$pVGRF_N_Neugebauer[i] <-
exp(5.247 + (0.271 * whole_sample_df$pVACC_g[i]) + (0.014 * whole_sample_df$body_mass[i]))
}Bland-Altman plots were used to confront the agreement between pVGRF measured with FP and those predicted using both regression equations.
Code to generate these plots can be seen here.
A lower dispersion around the bias value, as well as a lower percentage of values falling off the LoA can be appreciated in our equation. Additionally, while the bias in our equation tended to be zero, in the reference equation bias was always > 0 (p < 0.05), showing a consistent underestimation of pVGRF.
Also, to assess the prediction accuracy, MAE, MAPE and RMSE were calculated using pVGRF values predicted from both equations.
non_obese_our_accuracy <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_predicted")
non_obese_Neug_accuracy <- accuracy_indices(non_obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
obese_our_accuracy <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_predicted")
obese_Neug_accuracy <- accuracy_indices(obese_df, "pVGRF_N", "pVGRF_N_Neugebauer")
whole_sample_our_accuracy <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_predicted")
whole_sample_Neug_accuracy <- accuracy_indices(whole_sample_df, "pVGRF_N", "pVGRF_N_Neugebauer")The accuracy indices from our equation were all substantially lower compared to the reference equation indices, with an overall MAPE approximately 3 times smaller for our equation. For both equations, the MAPE was lower in the obese subsample than in the non-obese subsample
For the LR prediction, the same steps were followed as the GRF prediction. The code to generate the models can be seen bellow.
# For resultant peak loading rate
back_res_LR_LMM <- lme(
fixed = pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back_LR,
na.action = na.omit
)
r2_back_res_LR_LMM <- rsquared(back_res_LR_LMM)
summary(back_res_LR_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back_LR
## AIC BIC logLik
## 4874.519 4903.336 -2429.26
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.4054883 2372.028
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.6605835
## Fixed effects: pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -287.0209 1590.2580 211 -0.180487 0.8569
## pRATR_gs 572.7967 102.1134 211 5.609416 0.0000
## I(pRATR_gs^2) -9.8958 1.6666 211 -5.937732 0.0000
## body_mass 18.1178 18.8758 55 0.959841 0.3413
## pRATR_gs:body_mass 3.4078 1.1110 211 3.067340 0.0024
## Correlation:
## (Intr) pRATR_g I(RATR bdy_ms
## pRATR_gs -0.720
## I(pRATR_gs^2) 0.075 -0.336
## body_mass -0.950 0.603 0.146
## pRATR_gs:body_mass 0.693 -0.810 -0.243 -0.738
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.4916791 -0.5669550 -0.1458794 0.4684454 4.2366891
##
## Number of Observations: 271
## Number of Groups: 57
## Response family link method Marginal Conditional
## 1 pRLR_Ns gaussian identity none 0.7500468 0.7500468
hip_res_LR_LMM <- lme(
fixed = pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip_LR,
na.action = na.omit
)
r2_hip_res_LR_LMM <- rsquared(hip_res_LR_LMM)
summary(hip_res_LR_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip_LR
## AIC BIC logLik
## 4917.796 4946.584 -2450.898
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.5023566 2529.466
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.6017064
## Fixed effects: pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -3510.410 1846.9416 210 -1.900661 0.0587
## pRATR_gs 514.898 110.2864 210 4.668734 0.0000
## I(pRATR_gs^2) -8.639 1.4312 210 -6.035804 0.0000
## body_mass 51.937 20.6749 55 2.512072 0.0150
## pRATR_gs:body_mass 2.929 1.0337 210 2.833284 0.0051
## Correlation:
## (Intr) pRATR_g I(RATR bdy_ms
## pRATR_gs -0.785
## I(pRATR_gs^2) 0.329 -0.615
## body_mass -0.948 0.657 -0.096
## pRATR_gs:body_mass 0.748 -0.829 0.104 -0.780
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.22341333 -0.61951158 -0.07192759 0.43894165 3.05373725
##
## Number of Observations: 270
## Number of Groups: 57
## Response family link method Marginal Conditional
## 1 pRLR_Ns gaussian identity none 0.7021441 0.7021441
# For vertical peak loading rate
back_vert_LR_LMM <- lme(
fixed = pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = back_LR,
na.action = na.omit
)
r2_back_vert_LR_LMM <- rsquared(back_vert_LR_LMM)
summary(back_vert_LR_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: back_LR
## AIC BIC logLik
## 4858.839 4887.626 -2421.419
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.4732634 2393.903
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.6661911
## Fixed effects: pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -324.0761 1622.7623 210 -0.199706 0.8419
## pVATR_gs 552.8242 103.4180 210 5.345534 0.0000
## I(pVATR_gs^2) -11.9453 1.8892 210 -6.322881 0.0000
## body_mass 18.1405 19.3620 55 0.936915 0.3529
## pVATR_gs:body_mass 3.9586 1.1637 210 3.401656 0.0008
## Correlation:
## (Intr) pVATR_g I(VATR bdy_ms
## pVATR_gs -0.722
## I(pVATR_gs^2) 0.031 -0.287
## body_mass -0.949 0.604 0.193
## pVATR_gs:body_mass 0.689 -0.790 -0.326 -0.743
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.07632135 -0.62184789 -0.08981213 0.46058800 3.09934690
##
## Number of Observations: 270
## Number of Groups: 57
## Response family link method Marginal Conditional
## 1 pVLR_Ns gaussian identity none 0.7283041 0.7283041
hip_vert_LR_LMM <- lme(
fixed = pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass,
random = ~ 1 | ID,
method = "ML",
correlation = corAR1(),
data = hip_LR,
na.action = na.omit
)
r2_hip_vert_LR_LMM <- rsquared(hip_vert_LR_LMM)
summary(hip_vert_LR_LMM)## Linear mixed-effects model fit by maximum likelihood
## Data: hip_LR
## AIC BIC logLik
## 4907.858 4936.645 -2445.929
##
## Random effects:
## Formula: ~1 | ID
## (Intercept) Residual
## StdDev: 0.4267048 2637.244
##
## Correlation Structure: AR(1)
## Formula: ~1 | ID
## Parameter estimate(s):
## Phi
## 0.6724967
## Fixed effects: pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass
## Value Std.Error DF t-value p-value
## (Intercept) -2687.8662 1930.3344 210 -1.392436 0.1653
## pVATR_gs 407.8434 112.4430 210 3.627111 0.0004
## I(pVATR_gs^2) -7.6603 1.4186 210 -5.400053 0.0000
## body_mass 45.8905 21.7542 55 2.109503 0.0395
## pVATR_gs:body_mass 3.8995 1.0739 210 3.631057 0.0004
## Correlation:
## (Intr) pVATR_g I(VATR bdy_ms
## pVATR_gs -0.763
## I(pVATR_gs^2) 0.305 -0.596
## body_mass -0.951 0.643 -0.086
## pVATR_gs:body_mass 0.729 -0.839 0.101 -0.756
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.4906751 -0.5800428 -0.1047851 0.4384381 2.9373647
##
## Number of Observations: 270
## Number of Groups: 57
## Response family link method Marginal Conditional
## 1 pVLR_Ns gaussian identity none 0.6868396 0.6868396
In all models, conditional R2 values ranged from 0.69 to 0.75.
Similar to the GRF prediction models, LR model validation was assessed by the LOOCV method with the cross_validate_mixed_model() function and is shown below.
# For resultant peak loading rate
fix_eff <- pRLR_Ns ~ pRATR_gs + I(pRATR_gs^2) + body_mass + pRATR_gs:body_mass
rand_eff <- ~ 1 | ID
# Back
LOOCV_back_res_LR_LMM <- do.call(rbind, (lapply(unique(back_LR$ID), cross_validate_mixed_model, df = back_LR)))
# Hip
LOOCV_hip_res_LR_LMM <- do.call(rbind, (lapply(unique(hip_LR$ID), cross_validate_mixed_model, df = hip_LR)))
# For vertical peak loading rate
fix_eff <- pVLR_Ns ~ pVATR_gs + I(pVATR_gs^2) + body_mass + pVATR_gs:body_mass
rand_eff <- ~ 1 | ID
# Back
LOOCV_back_vert_LR_LMM <- do.call(rbind, (lapply(unique(back_LR$ID), cross_validate_mixed_model, df = back_LR)))
# Hip
LOOCV_hip_vert_LR_LMM <- do.call(rbind, (lapply(unique(hip_LR$ID), cross_validate_mixed_model, df = hip_LR)))Bland-Altman plots were also used to examine the agreement between actual and predicted pLR for both resultant and its vertical component.
Bland-Altman plots can be seen below. Panel C shows plots of pRLR and panel D shows plots of pVLR.
source(here("figs", "fig2.R"))
BA_plot_grid_1 <- plot_grid(
back_res_LR_BA_plot + theme(legend.position = "none"),
back_vert_LR_BA_plot + theme(legend.position = "none"),
hip_res_LR_BA_plot + theme(legend.position = "none"),
hip_vert_LR_BA_plot + theme(legend.position = "none"),
labels = c("C", "D", "", ""),
align = "h", vjust = 1, label_size = 16,
ncol = 2, nrow = 2
)
legend <- get_legend(back_res_LR_BA_plot)
BA_plot_grid <- plot_grid(BA_plot_grid_1, legend, ncol = 1, rel_heights = c(1, 0.1))
BA_plot_gridCode to generate these plots can be seen here.
One-sample t-tests were run to check whether bias from each accelerometer placement, in both pRGRF and pVGRF, was significantly different than 0. These tests were run using the t.test() function of the base R with the argument mu = 0. No significant differences were found.
##
## One Sample t-test
##
## data: LOOCV_back_res_LR_LMM$diff
## t = -0.95034, df = 269, p-value = 0.3428
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -418.0833 145.8665
## sample estimates:
## mean of x
## -136.1084
##
## One Sample t-test
##
## data: LOOCV_hip_res_LR_LMM$diff
## t = -0.73462, df = 269, p-value = 0.4632
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -418.3023 190.9683
## sample estimates:
## mean of x
## -113.667
##
## One Sample t-test
##
## data: LOOCV_back_vert_LR_LMM$diff
## t = -0.97483, df = 269, p-value = 0.3305
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -428.7289 144.7705
## sample estimates:
## mean of x
## -141.9792
##
## One Sample t-test
##
## data: LOOCV_hip_vert_LR_LMM$diff
## t = -1.1318, df = 269, p-value = 0.2587
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -497.1642 134.2107
## sample estimates:
## mean of x
## -181.4768
Also, linear regressions were applied to identify if there was any proportional bias, that is, if bias was related with the magnitude of the mean between measured and predicted pGRF (10). Linear regressions were run using the lm() function of the base R.
Results showed a significant proportional bias (p < 0.05), however, with a low magnitude. These results suggest that despite there is a trend for underestimation at increasingly higher pLR values, the magnitude of this effect is neglectable.
# For resultant peak loading rate
# Back
back_res_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_res_LR_LMM)
# Hip
hip_res_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_res_LR_LMM)
# For vertical peak loading rate
# Back
back_vert_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_back_vert_LR_LMM)
# Hip
hip_vert_LR_BA_plot_LR <- lm(diff ~ mean, data = LOOCV_hip_vert_LR_LMM)To evaluate models prediction accuracy, MAE, MAPE and RMSE were computed. These indices were computed with the accuracy_indices() function.
# For resultant peak loading rate
# Back
back_res_LR_accuracy <- accuracy_indices(LOOCV_back_res_LR_LMM, "pRLR_Ns", "pRLR_Ns_predicted")
# Hip
hip_res_LR_accuracy <- accuracy_indices(LOOCV_hip_res_LR_LMM, "pRLR_Ns", "pRLR_Ns_predicted")
# For vertical peak loading rate
# Back
back_vert_LR_accuracy <- accuracy_indices(LOOCV_back_vert_LR_LMM, "pVLR_Ns", "pVLR_Ns_predicted")
# Back
hip_vert_LR_accuracy <- accuracy_indices(LOOCV_hip_vert_LR_LMM, "pVLR_Ns", "pVLR_Ns_predicted")Table below shows these indices for each accelerometer placement in both pRLR and pVLR.
accuracy <- data.frame(
acc_placement = c(rep("Back", 2), rep("Hip", 2)),
condition = rep(c("pRLR", "pVLR"), 2),
MAE = c(
round(back_res_LR_accuracy[[5]], 1), round(back_vert_LR_accuracy[[5]], 1),
round(hip_res_LR_accuracy[[5]], 1), round(hip_vert_LR_accuracy[[5]], 1)
),
MAPE = c(
paste(round(back_res_LR_accuracy[[12]], 1), "%", sep = ""), paste(round(back_vert_LR_accuracy[[12]], 1), "%", sep = ""),
paste(round(hip_res_LR_accuracy[[12]], 1), "%", sep = ""), paste(round(hip_vert_LR_accuracy[[12]], 1), "%", sep = "")
),
RMSE = c(
round(back_res_LR_accuracy[[7]], 1), round(back_vert_LR_accuracy[[7]], 1),
round(hip_res_LR_accuracy[[7]], 1), round(hip_vert_LR_accuracy[[7]], 1)
)
)
kable(accuracy, col.names = c("Accelerometer Placement", "LR", "MAE", "MAPE", "RMSE"), align = c("l", "l", "r", "r", "r"))| Accelerometer Placement | LR | MAE | MAPE | RMSE |
|---|---|---|---|---|
| Back | pRLR | 1732.7 | 19.5% | 2352.9 |
| Back | pVLR | 1771.9 | 19.8% | 2393.0 |
| Hip | pRLR | 1866.8 | 20.6% | 2540.3 |
| Hip | pVLR | 1975.5 | 22.4% | 2636.1 |
Also, to assess if pLR predicted values significantly differed from the actual values measured by the FP, repeated measures ANOVA were run. As with the ANOVA for the pGRF, here walking speeds, accelerometer placements (lower back and hip), and the interaction effect (speed X accelerometer placements) were considered for analysis.
Code to tidy the data is shown below.
## Predicted pRLR
back_res_LR_pred <- LOOCV_back_res_LR %>%
select(ID, speed, pRLR_Ns_predicted) %>%
spread(key = speed, value = pRLR_Ns_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = back
)
hip_res_LR_pred <- LOOCV_hip_res_LR %>%
select(ID, speed, pRLR_Ns_predicted) %>%
spread(key = speed, value = pRLR_Ns_predicted) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = hip
)
res_pred_LR_df <- back_res_LR_pred %>%
full_join(hip_res_LR_pred, by = c("ID", "speed")) %>%
na.omit()
## Actual pRLR
back_res_actual_LR <- LOOCV_back_res_LR %>%
select(ID, speed, pRLR_Ns) %>%
spread(key = speed, value = pRLR_Ns) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_back
)
hip_res_actual_LR <- LOOCV_hip_res_LR %>%
select(ID, speed, pRLR_Ns) %>%
spread(key = speed, value = pRLR_Ns) %>%
na.omit() %>%
gather(
`2`, `3`, `4`, `5`, `6`,
key = speed,
value = actual_hip
)
res_actual_LR_df <- back_res_actual_LR %>%
full_join(hip_res_actual_LR, by = c("ID", "speed")) %>%
na.omit()
res_actual_LR_df <- res_actual_LR_df %>%
mutate(actual = (actual_back + actual_hip) / 2) %>%
select(ID, speed, actual)
## Merge predicted and actual data frames
res_ANOVA_LR_df <- res_pred_LR_df %>%
full_join(res_actual_LR_df, by = c("ID", "speed")) %>%
gather(
back, hip, actual,
key = "group",
value = "pRLR"
)
res_ANOVA_LR_df$ID <- as.factor(res_ANOVA_LR_df$ID)
res_ANOVA_LR_df$speed <- as.factor(res_ANOVA_LR_df$speed)
res_ANOVA_LR_df$group <- as.factor(res_ANOVA_LR_df$group)Repeated measures ANOVAs were performed using the ezANOVA() function of the ez package. Code to run the analysis is shown below.
##
## Pairwise comparisons using paired t tests
##
## data: s2_res$pRLR and s2_res$group
##
## actual back
## back 0.00016 -
## hip 0.00011 0.72304
##
## P value adjustment method: holm
These results can be seen in figure below.
## ─ Session info ──────────────────────────────────────────────────────────
## setting value
## version R version 3.5.2 (2018-12-20)
## os macOS 10.15.1
## system x86_64, darwin15.6.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/Lisbon
## date 2019-12-23
##
## ─ Packages ──────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 3.5.0)
## acepack 1.4.1 2016-10-29 [1] CRAN (R 3.5.0)
## assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.0)
## backports 1.1.2 2017-12-13 [1] CRAN (R 3.5.0)
## base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.5.0)
## broman * 0.68-2 2018-07-25 [1] CRAN (R 3.5.0)
## broom 0.4.5 2018-07-03 [1] CRAN (R 3.5.0)
## callr 2.0.4 2018-05-15 [1] CRAN (R 3.5.0)
## car 3.0-2 2018-08-23 [1] CRAN (R 3.5.0)
## carData 3.0-2 2018-09-30 [1] CRAN (R 3.5.0)
## cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.5.0)
## checkmate 1.8.5 2017-10-24 [1] CRAN (R 3.5.0)
## cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.0)
## cluster 2.0.7-1 2018-04-13 [1] CRAN (R 3.5.2)
## colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.0)
## cowplot * 0.9.3 2018-07-15 [1] CRAN (R 3.5.0)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.0)
## curl 3.2 2018-03-28 [1] CRAN (R 3.5.0)
## data.table 1.12.0 2019-01-13 [1] CRAN (R 3.5.2)
## desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.0)
## devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.0)
## digest 0.6.19 2019-05-20 [1] CRAN (R 3.5.2)
## dplyr * 0.8.0.1 2019-02-15 [1] CRAN (R 3.5.2)
## evaluate 0.10.1 2017-06-24 [1] CRAN (R 3.5.0)
## ez * 4.4-0 2016-11-02 [1] CRAN (R 3.5.0)
## fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.0)
## forcats * 0.3.0 2018-02-19 [1] CRAN (R 3.5.0)
## foreign 0.8-71 2018-07-20 [1] CRAN (R 3.5.2)
## Formula 1.2-3 2018-05-03 [1] CRAN (R 3.5.0)
## fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.0)
## ggplot2 * 3.1.0 2018-10-25 [1] CRAN (R 3.5.0)
## glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.0)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 3.5.0)
## gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.0)
## haven 2.0.0 2018-11-22 [1] CRAN (R 3.5.0)
## here * 0.1 2017-05-28 [1] CRAN (R 3.5.0)
## highr 0.7 2018-06-09 [1] CRAN (R 3.5.0)
## Hmisc 4.3-0 2019-11-07 [1] CRAN (R 3.5.2)
## hms 0.4.2 2018-03-10 [1] CRAN (R 3.5.0)
## htmlTable 1.12 2018-05-26 [1] CRAN (R 3.5.0)
## htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.0)
## htmlwidgets 1.3 2018-09-30 [1] CRAN (R 3.5.0)
## httr 1.3.1 2017-08-20 [1] CRAN (R 3.5.0)
## jsonlite 1.5 2017-06-01 [1] CRAN (R 3.5.0)
## knitr * 1.21 2018-12-10 [1] CRAN (R 3.5.0)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
## lattice 0.20-38 2018-11-04 [1] CRAN (R 3.5.2)
## latticeExtra 0.6-28 2016-02-09 [1] CRAN (R 3.5.0)
## lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.0)
## lme4 1.1-19 2018-11-10 [1] CRAN (R 3.5.0)
## lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.5.0)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.0)
## MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.2)
## Matrix 1.2-15 2018-11-01 [1] CRAN (R 3.5.2)
## memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.0)
## mgcv 1.8-26 2018-11-21 [1] CRAN (R 3.5.2)
## minqa 1.2.4 2014-10-09 [1] CRAN (R 3.5.0)
## mnormt 1.5-5 2016-10-15 [1] CRAN (R 3.5.0)
## modelr 0.1.2 2018-05-11 [1] CRAN (R 3.5.0)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.0)
## nlme * 3.1-137 2018-04-07 [1] CRAN (R 3.5.2)
## nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.5.0)
## nnet 7.3-12 2016-02-02 [1] CRAN (R 3.5.0)
## openxlsx 4.1.0.1 2019-05-28 [1] CRAN (R 3.5.2)
## pbkrtest 0.4-7 2017-03-15 [1] CRAN (R 3.5.0)
## piecewiseSEM * 2.0.2 2018-07-24 [1] CRAN (R 3.5.0)
## pillar 1.3.1 2018-12-15 [1] CRAN (R 3.5.0)
## pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.0)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.0)
## pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.0)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.0)
## prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.0)
## processx 3.1.0 2018-05-15 [1] CRAN (R 3.5.0)
## psych 1.8.4 2018-05-06 [1] CRAN (R 3.5.0)
## purrr * 0.3.0 2019-01-27 [1] CRAN (R 3.5.2)
## R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.2)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.5.0)
## Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.0)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.5.0)
## readxl 1.1.0 2018-04-20 [1] CRAN (R 3.5.0)
## remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.0)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.0)
## rio 0.5.16 2018-11-26 [1] CRAN (R 3.5.0)
## rlang 0.4.1 2019-10-24 [1] CRAN (R 3.5.2)
## rmarkdown 1.10 2018-06-11 [1] CRAN (R 3.5.0)
## rpart 4.1-13 2018-02-23 [1] CRAN (R 3.5.2)
## rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.0)
## rstudioapi 0.7 2017-09-07 [1] CRAN (R 3.5.0)
## rvest 0.3.2 2016-06-17 [1] CRAN (R 3.5.0)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.0)
## sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.0)
## stringi 1.3.1 2019-02-13 [1] CRAN (R 3.5.2)
## stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.5.2)
## survival 2.43-3 2018-11-26 [1] CRAN (R 3.5.2)
## testthat 2.0.0 2017-12-13 [1] CRAN (R 3.5.0)
## tibble * 2.0.1 2019-01-12 [1] CRAN (R 3.5.2)
## tidyr * 0.8.2 2018-10-28 [1] CRAN (R 3.5.0)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.0)
## tidyverse * 1.2.1 2017-11-14 [1] CRAN (R 3.5.0)
## usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.0)
## utf8 1.1.4 2018-05-24 [1] CRAN (R 3.5.0)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.0)
## xfun 0.4 2018-10-23 [1] CRAN (R 3.5.0)
## xml2 1.2.0 2018-01-24 [1] CRAN (R 3.5.0)
## yaml 2.1.19 2018-05-01 [1] CRAN (R 3.5.0)
## zip 1.0.0 2017-04-25 [1] CRAN (R 3.5.0)
##
## [1] /Library/Frameworks/R.framework/Versions/3.5/Resources/library
1. Guthold R, Stevens GA, Riley LM, Bull FC. Worldwide trends in insufficient physical activity from 2001 to 2016: a pooled analysis of 358 population-based surveys with 19 million participants. Lancet Glob Heal. 2018;6(10):e1077–86.
2. Li X, Gong X, Jiang W. Abdominal obesity and risk of hip fracture: a meta-analysis of prospective studies. Osteoporos Int. 2017;28(10):2747–57.
3. Hind K, Rudman H, Pearce M, Treadgold L, Birrell F. Obesity, bone density for weight and prevalent vertebral fracture at age 62 years: the Newcastle Thousand Families Study. J Clin Densitom. 2018;21(4):606.
4. Kaze AD, Rosen HN, Paik JM. A meta-analysis of the association between body mass index and risk of vertebral fracture. Osteoporos Int. 2018;29(1):31–9.
5. Viester L, Verhagen EA, Hengel KM, Koppes LL, Van Der Beek AJ, Bongers PM. The relation between body mass index and musculoskeletal symptoms in the working population. BMC Musculoskelet Disord. 2013;14(238).
6. Nakagawa S, Schielzeth H, O’Hara RB. A general and simple method for obtainingR2from generalized linear mixed-effects models. Methods in Ecology and Evolution. 2013;4(2):133–42.
7. Jaeger BC, Edwards LJ, Das K, Sen PK. An r2 statistic for fixed effects in the generalized linear mixed model. Journal of Applied Statistics. 2016;44(6):1086–105.
8. Staudenmayer J, Zhu W, Catellier DJ. Statistical considerations in the analysis of accelerometry-based activity monitor data. Medicine & Science in Sports & Exercise. 2012;44(1 Suppl 1):S61–7.
9. Bland JM, Altman DG. Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet. 1986;1(8476):301–10.
10. Giavarina D. Understanding bland altman analysis. Biochemia Medica. 2015;25(2):141–51.
11. Neugebauer JM, Collins KH, Hawkins DA. Ground reaction force estimates from ActiGraph GT3X+ hip accelerations. PLoS One. 2014/06/11. 2014;9(6):e99023.